function B_rho=B_coil_rho(z_,rho_,I_,a_)

   
    %constant
   mu_0=1.2566370614e-6;
        
    k=sqrt(4*a_.*rho_./((rho_+a_).^2+z_.^2));
    k2=k.^2;
    [K,E] = ellipke(k2); %calculates the elliptic integrals of first and second kind
    
    aux1=mu_0*I_/(2*pi);
    aux2=1./sqrt((rho_+a_).^2+z_.^2);
    aux3=1./(a_-rho_).^2+z_.^2;
    aux4=mu_0*I_*(a_.^2)/2;
        
    B_rho=aux1.*(z_./rho_).*aux2.*((a_^2+rho_.^2+z_.^2).*aux3.*E-K);
    B_z=aux1.*aux2.*((a_.^2-rho_.^2-z_.^2).*aux3.*E-K);
    
    %replacing NaNs anf Infs
    %
    B_z(isnan(B_rho))=aux4.*(z_(isnan(B_rho)).^2+a_^2).^(-3/2);
    B_z(isinf(B_rho))=aux4.*(z_(isinf(B_rho)).^2+a_^2).^(-3/2);
    B_rho(isnan(B_rho))=0;
    B_rho(isinf(B_rho))=0;
   %
    B_=B_rho;
end